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ABSTRACT 


A numerical  method  for  solving  two-dimensional  boundary- 
value  problems  related  to  time-harmonic  potential  flows  with 
a free-surface  is  introduced  in  this  paper.  There  is  no 
restrictions  on  the  body  geometry  nor  the  bottom  topography. 
The  entire  fluid  domain  is  subdivided  into  two  regions,  an 
inner  region  where  all  the  geometrical  changes  occur  and 
an  outer  region  where  the  fluid  depths  are  constant.  Appli- 
cation of  Green's  theorem  for  an  inner  region  results  in  an 
integral  relation  while  the  outer  region  is  described  by 
eigenfunction  expansions.  A continuity  of  pressure  and  nor- 
mal velocity  across  junction  boundary  imposes  an  unique  solu- 
tion in  the  entire  fluid  domain.  The  method  is  applied  to 
both  radiation  and  scattering  problems  and  test  results  in 
water  of  finite  or  infinite  depth  agree  well  with  the  exis- 
ting results. 
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I 


INTRODUCTION 


I . 


Hydrodynamic  research  of  two-dimensional  bodies  har- 
monically oscillating  in  or  below  the  free  surface  of  an 
ideal  fluid  has  increased  in  importance  in  past  decades  and 
has  been  studied  by  a number  of  researchers.  The  'radiation' 
problem  applies  to  situation  when  waves  are  generated  by  for- 
ced oscillatory  motion  of  the  body  in  an  otherwise  undis- 
turbed water.  The  modern  history  of  this  subject  began  with 
Ursell  (1949) , who  formulated  and  solved  the  boundary-value 
problem  for  a semi-immersed  heaving  circular  cylinder  using 
linearized  free-surface  theory.  He  represented  the  velocity 
potential  as  the  sum  of  an  infinite  set  of  multipoles,  each 
satisfying  the  linear  free-surface  condition  and  each  being 
multiplied  by  a coefficient  determined  by  requiring  the 
series  to  satisfy  the  kinematic  boundary  condition  at  a 
number  of  points  on  the  cylinder.  Grim (1953)  used  a varia- 
tion of  the  Ursell  method  to  solve  the  problem  for  two- 
parameter  Lewis-form  cylinders  by  conformal  mapping  onto 
circle.  Tasai(1959)  and  Porter (1960) , using  the  Ursell 
approach,  obtained  the  added  mass  and  damping  for  oscillating 
contours  mappable  onto  a circle  by  the  more  general  Theodorsen 
transformation.  Ogilvie (1963)  calculated  the  hydrodynamic 
forces  on  completely  submerged  circular  cylinder.  Lebreton 


8 


and  Margnac (1966) , Frank  (1967)  distributed  the  so-called 
wave-sources  on  the  body  contour  and  obtained  the  strength 
of  these  sources  by  solving  a Fredholm  integral  equation  of 
the  second  kind.  This  approach  is  still  applicable  when 
the  depth  is  net  constant,  in  which  case  the  unknown  strength 
in  the  integral  equation  will  include  also  contribution  from 
the  bottom  boundary. 

The  'scattering'  problem  customarily  refers  to  the 
interaction  between  incident  plane  waves  and  a fixed  body. 
Various  approaches  have  been  used  to  tackle  these  problems 
associated  with  simple  geometries.  Newman (1965)  provided 
results  for  the  reflection  and  transmission  coefficients 
over  an  infinite  step  by  matching  wave-maker  solutions. 

Miles  (1967)  showed  a more  exact  solution  of  the  unequal- 
depths  problem  by  using  Schwinger's  variational  formulation. 
Hilaly(1967)  treated  the  same  problem  by  matching  eigenfunc- 
tion expansions.  A more  elegant  approach,  based  upon  a 
single  source  solution  valid  everywhere,  may  be  constructed 
from  the  Green  function  derived  by  Evans (1971).  Recently, 
Bai(1972)  and  Yeung  (1973)  investigated  two  new  approaches 
for  treating  the  problems  of  general  bottom  topography.  The 
former  utilized  a variational  formulation  with  the  fluid 
domain  represented  by  finite  elements.  The  latter  approach, 
on  the  other  hand,  consisted  of  applying  a simple  source 
function  to  the  fluid  domain,  resulting  in  an  integral 


equation  along  the  domain  boundary  that  has  to  be  solved. 

While  from  the  physical  viewpoint,  the  scattering  and 
radiation  problems  appear  to  be  unrelated,  from  the  mathma- 
tical  viewpoint  they  differ  only  in  terms  of  the  boundary 
condition  on  the  body  and  the  far-field  conditions.  Many 
useful  relations  between  them  can  be  obtained  as  a results 
of  application  of  Green's  Theorem  or  a consequence  of 
reciprocity  principle. 

The  method  investigated  in  this  thesis  stems  from  Yeung 
(1973)  and  Bai  and  Yeung  (1974).  The  essence  of  the  present 
approach  lies  in  the  exploitation  of  more  than  a single  term 
in  an  eigenfunction  expansions  as  previously  used  by  Yeung 
(1973).  TniL  will  have  the  net  effect  of  shortening  the 
boundary  of  the  fluid,  hence  leads  to  an  increase  in  compu- 
tational efficiency.  The  overall  approach  to  the  solution 
of  the  two-dimensional  radiation  or  scattering  problems  may 
be  summarized  as  follows.  The  entire  fluid  domain  is  sub- 
divided into  two  regions,  an  inner  region  where  all  the 
geometrical  changes  occur  and  an  outer  region  (s)  where  the 
fluid  depth  is  a constant.  For  the  inner  region,  application 
of  Green's  second  identity  to  the  unknown  velocity  potential, 
Cp(p)  and  the  simple  source  function  log(l/r)  give  us  an 
integral  relation  between  cp (p ) and  its  boundary  values. 

For  the  outer  region,  the  velocity  potential  is  written  in 
the  form  of  an  eigenfunction  expansion  with  unknown  complex 
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coef f icients . At  the  junction  boundaries  where  the  inner 
and  outer  regions  meet,  a matching  of  the  velocity  potential 
and  normal  velocities  is  performed.  This  implies  physically 
a continuity  of  pressure  and  flux  across  the  artificial 
junction  or  "radiation"  boundary.  Such  a matching  process 
determines  the  solution  uniquely  in  the  entire  fluid  domain. 
We  found  that  this  method  not  only  gives  very  accurate 
results  compared  with  known  values  but  also  reduces  compu- 
tational time  remarkablely . 
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II.  FORMULATION  OF  THE  PROBLEM 

Consider  the  two-dimensional  time-harmonic  motion  of 
an  ideal  fluid  with  a f ree-surf ace , as  shown  in  Figure  1 on 
page  17.  A right-handed  rectangular  coordinate  system  is 
used.  The  x-axis  coincides  with  the  undisturbed  f ree- 
surf  ace,  the  positive  y-axis  is  taken  upwards.  The  depth 
of  fluid  can  be  either  finite  or  infinite.  In  the  finite 
case,  the  bottom  topography  can  be  arbitrary,  but  must 
approach  asymptotically  constant  depths  h+  and  h , which 
need  not  be  equal.  In  order  to  use  linearized  water-wave 
theory,  the  following  assumptions  are  required: 

a) ,  the  fluid  is  inviscid  and  incompressible; 

b)  . the  effects  of  surface  tensions  are  negligible; 

c)  . the  fluid  flow  is  irrotational ; 

d)  . the  fluid  motion  amplitudes  and  velocities  are 
small  enough  that  all  but  linear  terms  of  the  free-surface 
condition  and  the  Bernoulli  equation  are  neglected.  The 
body  boundary  condition  can  be  satisfied  at  an  "mean"  posi- 
tion, if  the  oscillations  performed  by  the  body  are  assumed 
small.  Under  the  following  assumptions,  there  exist  a 
time-harnomic  velocity  potential  which  describes  the  motion 
of  fluid. 


II. 1 Governing  Equations  and  Boundary  Conditions 


Let  §i.(x,y,t)  be  the  velocity  potential  describing 
the  flow  field.  The  continuity  equation  requires  ^ to 
satisfy  Laplace  equation.  Let  cr  be  the  angular  frequency 
of  the  time-harmonic  solution,  then  we  have: 

§>(x,^,t)  = Re  ( jp  d(t)]  (2.1) 

where  cp  = CP, is  the  unknown  complex-valued  spatial 
function,  and  CL( t)  , the  complex  time-harmonic  amplitude  of 
body,  i.e. 

CX  ( Jr)  = Re  [ ftU)l  = Re  [ ($i  + < $2 ")  , si- JT  (2.2) 

where  cX (t)  is  the  body  motion  amplitude. 

The  unknown  complex-valued  spatial  function,  £j?(x,y)  , 
must  satisfy  Laplace  equation  throughout  the  fluid  domain; 


( /x,  jp  = o 


(2.3a) 


where  v = 


2 _ 


e-*1 


+ 


* r)<J> 

= a* 

= if 


and  also  the  following  linerized  boundary  conditions; 

= 0 on  SF  (2.3b) 

- 0 on  S>B  (2.3c) 

(Pn  = Vn  = f C S)  2>0  (2.3d) 

where  y = -li—  , \ being  the  acceleration  of  gravity. 
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Here,  Sp  is  the  undisturbed  free  surface,  the  bottom 
surface,  and  SQ  the  wetted  surface  of  the  body  at  its 
equilibrium  position.  The  normal  vector  rf  is  taken  to  be 
positive  when  pointing  out  of  the  fluid  domain  with  compon- 
ents (n  ,n  ).  The  only  non-homogeneous  boundary  condition 
x y 

is  (2.3d),  which  is  given  by  f(s)=n  , n , and  (xn  -yn  ), 

x y y x 

each  corresponding  to  the  sway,  heave  and  roll  motion  of 
the  body  respectively.  The  radiation  condition,  which  states 
that  the  disturbance  sufficiently  far  from  the  generating 
points  must  be  outgoing  waves,  will  render  the  solution 
unique. 

11.2  Pressure,  Force,  Moment  and  Free-Surface  Elevation 

The  linearized  hydrodynamic  forces  and  moments  acting 
on  the  body  are  given  by; 

p = 1 fn  dA  (2.4a) 

Sb 

//|  = ( -p(' rx~rT)  d-4  (2.4b) 

k 

where  cU  is  the  infinitesimal  arc-length  element  and  "f7 
can  be  obtained  from  the  time -dependent  Bernoulli  equation, 

f " ffj  (2.4c) 

where  ^ denoting  the  density  of  fluid.  The  first  term 
in  (2.4c)  corresponds  to  the  hydrodynamic  pressure  which  is 
our  primary  concern  and  the  second  the  hydrostatic  pressure. 
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If  we  make  use  of  (2.2)  in  (2.4c),  (2.4a)  and  (2.4b)  then 


L 

a i 


yield , 


Fj  = “ ) U ~ 


-i-i 


nid-4  - fcrZ  <P3 

L i 4 11 


} ~ * » ^ 


(2.5a) 


M;  = -(  f^(Fxn)^  = -fi  -r^i.  [(  (p‘°n;^  , 

i V <5  ^lli.  <*  J s.  * 

■ - (2.5b) 

<f  "3 

where  i=l,2,3  correspond  to  the  sway,  heave  and  roll  motion 


respectively.  Let's  introduce  the  notation  F^M-j  and  define 

^ +<'ffV  = +.j=i,z, 3 (2.6) 

Se  <* 

then  the  force (or  moment)  component  in  the  j-th  direction 
can  be  written  as 

Fj  = (2,7) 

where  repeated  indices  are  summed  from  1 to  3.  The  are 

called  'added  mass'  or  'moment  of  added  mass',  depending  on 
their  dimensions,  and  each  represents  the  fluid  reaction 
in  the  j-th  direction  due  to  an  acceleration  of  the  body  in 
the  i-th  direction.  The  7\^  are  called  'damping  coeffici- 
ents' and  correspond  to  fluid  reactions  in  phase  with 
velocities.  The  free-surface  elevation  C,(x,t)  can  be 
obtained  from  the  dynamic  boundary  condition  of  the  free- 
surface  and  can  be  written  as 

- 15  - 
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III . 


SOLUTION  OF  THE  PROBLEM 


One  difficulty  in  the  numerical  solution  of  the  radia- 
tion or  forced-motion  problem  is  that  the  radiation  condition 
in  principle  should  be  applied  as  x=ioo  . In  practice,  the 
infinite  boundaries  have  to  be  truncated  at  a finite  range 
with  approprite  reasonings,  the  result  of  such  truncation 
implies  local  disturbances  must  be  accounted  for.  The 
behavior  of  local  disturbances  can  be  examined  by  making 
use  of  the  eigenfunction  expansions. 

First  let  Z+,Z_  be  two  arbitrarily  chosen  vertical 
boundaries  which  subdivide  the  fluid  domain  into  an  inner 
region  of  arbitrary  bottom  and  the  body  and  two  outer  regions 
of  horizontally  uniform  depths.  Let  these  regions  be  desig- 
nated I,  II  and  III  as  shown  in  Figure  1. 


! 
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Fig . 1 Coordinate  System  and  Cylinder  in  Oscillation 
- 17  - 


Let  and  % denote  the  complex-valued  potential 

functions  in  each  region,  then  we  have  the  boundary-value 


problems  as  follows: 


v‘cp 

- 0 

in 

fluid 

domain (I , 

II,  HI) 

V 

vcp  - o 

on 

SF 

(I, 

II, III) 

£ 

ii 

II 

£ 

on 

S 

o 

(I) 

ii 

o 

on 

SB 

(I, 

II, III) 

ju 

(<P*  + = o 

, "X  — — ± °° 

(II 

,111) 

We  need  some  additional  conditions  for  matching  the  pressure 
and  flux  which  are  proportional  to  the  velocity  potential 
and  normal  velocity  respectively,  at  the  arbitrarily  chosen 
boundaries; 

CPi  — tfz  7 <^>,n  ~ ~ 

(3.2) 

CP,  = % > <fU  = of 

With  these  matching  conditions,  the  complex-valued  functions 
($  ^ cpz  and  (Q  would  be  properly  coupled  together  and 
represent  a unique  solution  for  the  entire  fluid  domain. 


III.l  Application  of  Green's  Theorem 

Referring  to  Figure  1,  we  assume  that  Cp  is  a smooth 
function  whose  first  and  second  derivatives  are  continuous 
everywhere  in  the  fluid  domain.  Application  of  Green's 
second  identity  to  the  unknown  complex  function  and  the 
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source  function  log(l/r)  in  Region  I,  yields  the  following 
well-known  theorem  : 


2 TT  (p,  ( P)  = (3.3) 

where  P = (x,y)  is  a field  point  in  Region  I and  Q = (i-,']) 
is  a variable  point  on  the  fluid  boundary,  S=  souspiv2.+  usBiv2i_* 
Here  r is  the  distance  between  points  P and  Q,  ds (Q)  an 
infinitesimal  arc-length  element  along  the  fluid  boundary  S, 
and  ^ , fj  being  variables  of  integration  along  S. 

From  the  boundary  conditions  (3.1),  we  observe  that 
the  normal  derivatives  of  CR  are  either  known  or  expressed 
in  terms  of  ^ itself  on  S except  for  2+  and  2.  • Therefore 
(3.3)  can  be  written  as  follows  : 

ztt ( P)  = ^ cR(a)^ri^r  ci4(Gi) -f- | Cp,  ^Jojr-y^r  d& 

+ i2 1<P, i - <Pm  ^ <P,  ^ <3  • 4 > 

Note  that  we  have  place  no  restrictions  on  the  shape  of 
body  and  geometry  of  bottom.  This  integral  relation  merely 
represents  the  distributions  of  simple  source  and  normal 
dipoles  with  a certain  amount  of  strength  along  the  contour 
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III. 2 Eigenfunction  Expansions 

Now  let  us  consider  Regions  II  and  III.  The  boundary 
of  Region  II  consists  of  a free  surface  Sp2,  the  horizontal 
bottom  SB2  and  the  "junction"  boundary  2+  . By  using  the 
method  of  seperation  of  variables,  one  can  obtain  the  fun- 
damental solutions  in  the  form  of  eigenfunction  expansions 
easily. (Wehausen  & Laitione , 1960 ) 

These  eigenfunctions  are 

( co-ail  Wo (i+fl ) ^ coo 

| Co^JL  YY\t , £»<3  7>V R ) 

where  the  eigenvalues  mQ,  m^ , m2  , . . . , m^,...  are  roots  of 
the  transcendental  equations; 

W0  — 4 

t (3.5) 

inK  Aom.  ( K = I,  Z,  3,  • • - ■ ) 


Eigenvalues  are  shown  graphically  below  : 


The  orthogonality  of  eigenfunctions  can  be  easily  verified 
by  direct  computation  on  the  interval  -hiy  < 0.  These  are 
all  the  consequence  of  Sturm-Liouville  theory  (Hildebrand, 
1948).  With  a similar  representation  in  Region  III,  we 
may  write  the  solutions  as  follows; 

<P  (*,*/)  = coo 

1 " ° Coo^.  -ft."*  ic-i  ^ 

for  Regionll  (3.6) 


3;i:)  + z c;  e*’ 

CxmA  m0-h 


K-*-« 


doo  no« 


for  Region  111(3.7) 


where  the  superscripts  + and  - denote  Regions  II  and  III 
respectively.  The  coefficients  Aq  and  are  unknown 
complex  values.  Note  that  the  first  term  of  (3.6)  or  (3.7) 
corresponds  to  propagating  waves,  the  second  to  purely  local 
disturbance.  Note  also  the  functional  forms  of  (3.6)  and 
(3.7)  are  chosen  so  that  CP2  and  represent  outgoing 

waves  propagating  to  the  right  and  left  respectively.  It 
is  clearly  not  practical  to  take  an  infinite  number  of 
terms  in  (3.6)  and  (3.7).  To  assists  us  in  deciding  the 
number  of  terms  to  be  taken  in  the  series,  let  us  define  the 
decay  factor  from  (3.6)  or  (3.7)  by 

St*,*)  = e?’v< 


21 


Assuming  C^=  0(1),  and  that  E is  a tolerance  parameter,  say, 
”6  R L 

10  , we  may  choose  N ' , the  number  of  terms  to  be  taken 

in  the  series,  either  by  finding  NR,L  such  that  S (NR,  (jR)  <fc 

and  S (NL,^L)<  & for  a given  ^R  and  ^ ^ , o£  finding  ^R 

t R L 

and  satisfying  the  same  criteria  for  given  N ' . This 

serves  somewhat  as  a quide  only,  because  in  actuality  Ck 

could  be  quite  large.  However  convergence  of  the  series  is 

guaranteed  by  the  fact  that  C^e  mkx  remain  small. 

Now  we  are  in  a position  to  apply  the  matching  condition 

(3.2)  at  the  junction  boundaries  . if  we  make  use 

of  ££  and  % and  substitute  them  into  at  Z+  and  , 

we  may  rewrite  (3.4)  as  : 


where  the  functions  and  Giffr'^.J) 

by  the  following  integrals  ; 


are  defined 


(3.10) 


?K  j ^rVl  K , , 2 

GyP-VJ)  'k  ' t wJ°V  ) ’ 

From  now  on,  subscript  ' 1*  wil  be  dropped  with  the  under- 
standing that  we  seek  only  (f  in  Region  I. 

One  can  easily  recognize  that  the  last  four  terms 
of  (3.8)  represent  vertical  distributions  of  simple  sources 
and  normal  dipoles  with  the  functional  forms  of  the  strength 
being  the  same  as  the  eigenfunction  in  the  vertical  direction. 
Referring  to  (3.8)  we  see  that  the  spatial  velocity  potential 
^P(P)  will  be  known  everywhere  if  £j?(Q)  is  determined  on 

the  boundaries  S S„,  and  if  the  coefficients  of  the 

°u  Bu  E 

eigenfunctions  are  determined.  To  obtain  an  equation  for 
the  (P  ( Q ) , we  let  the  field  point  P approach  the  fluid 
boundary  S.  The  following  integral  equation  then  results  : 
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for  Pe  S. 


(3.11) 


Note  that  to  determine  the  unknown  complex  coefficients  of 
eigenfunction,  we  need  (NR+1)  distinct  locations  on  2+  and 
(N  +1)  distinct  locations  on  Z_  when  the  field  point  P is 
placed  at  each  junction  boundary.  In  this  case  CP(P)  can 
be  expressed  as  the  eigenfunction  expansions  themselves, i. e. , 


t L e 


+ j <P(|;A)r -vAji-)* 

, £>F 

_ CoA 


+ ZC*e 


_jSsML  + Cry.v  }-K^ 


Coo  vnK 


tzQe 

K-l 


So  „ , 


for  P € 2i+  (3.12) 


and 


SoySg 


t A e 
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> Ur  -t-  -m3  + 

\*  zc* e V- 


-A  e 


-3 


T +^r'(o,p 

f)  ~ 0 - ' o u 

rv)0  -K 


NlL  _ 

+ 2 Ck  e 


K“J 


_ tt  6^  rr?K  l&i2  +fZ(o,n)l  ==•  f -|(S)  Zt t~ 

coo  ryi-  -T  3 d l a 


for  P € 2)  (3.13) 
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The  method  of  discretization  will  be  used  to  obtain  an 
approximate  solution  of  the  integral  equation  (3.11)  for 
the  unknown  function  (P  ( Q)  . Such  method  reduces  the  integral 
equation  to  a set  of  linear  algebraic  equations.  We  follewed 
the  scheme  used  by  Yeung  (1973).  Referring  to  Figure  2 


Ng  : number  of  intervals  on  body 
N ?-N s'-  number  of  intervals  on  free-surfaoe 
eigenfunction-expansions  at  R.H.S. 

etc 


Fig.  2 Subdivision  of  Contour 


the  contour  S , S_and  S„  in  Region  I is  first  subdivided 
o F B 

into  a number  of  small  segments.  Then  each  segment  can 

be  considered  as  a straight  line  joining  a pair  of  grid 

points  ( ) , ( $.+l  , ) • 

For  a sufficiently  fine  division,  the  unknown  ty(Q) 

between  a pair  of  grid  points  may  be  approximated  by  some 

distribution  function.  In  present  case  <p(Q)  is  assumed  to 

be  a constant  function  and  represented  by  = cp ( ^V1-  » 

As  an  illustration,  let  us  consider  the  body  contour  integral 

along  S ; 
o 


I 


= Z . <?  TF 


where  SQis  broken  up  into  a number  of  small  segments  S^. 
in  each  segment,  (f(s)  may  be  approximated  by  constant. 
Hence  it  can  be  written  as  follows  ; 

S - f v-'j 

£>o  * 


with 


Qii 


IW 


Similar  operations  on  other  parts  of  the  contour  then  yield 


the  following  discretized  representation  of  integral  of  (3.11) 
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- TT  cp  ( X;  •>$.<: ) + 2,  Q[  Q^;  + Z Q:  (Qii~  wP-tf)  ■+•  A„  Se 

J-si  <r  <r  J ^ * 

Nr  -f  + ^4  - - 

+ 2 Q SK  + .2;  Q ;Q ^ + A0  S^k  + ^ Ck  SA.K  -+  Z.  ^ (Qtj  - ) 

K=t  1=  hL+l  * & K=l  * 6 * 


-i  ^c*i’W  ^‘J 


i=l , 2 , . . .Ng , Ng+1,...N2 


i=N3+l,  ...N4, 
i=N5+l,  ...N5, 


(3.14) 


where  (x^,y^)  denotes  the  control  or  observation  point  and 


IV’1 y 


p-j  ~ 


<ir1 i> 


x 

Qii-(  ik 

J K yV  > 


(3.15) 


S-o  = e 


(3.16a) 


± 

a,  = 


3k  ^ ~l^K  ^ ),  k=,’2-~ 


(3.16b) 


with  JK  and  (jy  defined  by  (3.9)  and  (3.10),  and  + ,- 
denoting  right  and  left  junction  or  "radiation"  boundary 
respectively.  Using  a more  compact  notation,  equation  (3.14) 


can  be  rewritten  as  follow 
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r 


V-  I 

4; 


Nf? 

% (-^%  +cGj  ) <Pj  + Al S7o  + 2 Ck  s7<  + Ao  Szo 

+ = i=l,2,  ...,N  (3.17) 

k>  I J 1 


where  is  the  Kronecker  delta  and  the  coefficient  G^j 

are  defined  as  follows 

Qn 


■j  or  ^ - I . ' ‘ ' • 5 ^34  h ' ‘ * ’ "> 

Qi!  - jor  jsKj-H,— .INI,  ; — . NT 


(3.18) 


When  we  place  the  control  points  on  junction  boundary  Z+  , 
(3.17)  needs  a slight  modification  for  (3.12)should  be  used 
instead  of  (3.11).  In  this  case,  the  unknown  velocity  poten- 
tial can  be  represented  by  eigenfunction  itself,  namely, 

a,rE)  - a*jsA^i± j:  — ^ 

cod  Cocm^-i 


for  P = ( $R,y)  (3.19) 


If  we  make  use  of  (3.19)  in  (3.12),  the  following  equation 
then  results: 


sis  .+  __■+  ^ /'-*■  -r  + 

Z <pt  % t X cpf  (Or;  - 1-  L Txo  + £ c»<  Uk 

J=l  i o j*^41  o * * 


:s;.  + xc;s = z %!  <*}•*,) 

for  i=N2+l , • • • ,N^  (3.20) 


4 Z cp.  <3  • • + A 
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where  are  defined  by 


To  = -f-  ^ &U°’  % > 

*°  m0+£+  ° o 

T- T = - TT  ***  "vjp^  + GtkCo’P 

c<sm  m?  -£+  K o 


(3.21) 


If  the  control  points  are  placed  on  Z>_  , the  velocity  poten- 
tial can  be  written  as  ; 


m0  A k'=' 


t&i  rt)*(p£~)  eW(<  T 
DIkT 


for  P = ( jj^y)  (3.22) 


and  the  following  results  are  obtained  by  making  use  of  (3.22) 
in  (3.13): 

Z 4>:  GU  + X CP-  ( Q • • - vP. . ) + A.  $10  + X 6*  S;k 

:=.!  <7  jcNlj-H  * 0 K ' 


+ Z<%<P,-  tltT,;  =2&}ft.3;) 

jrfij+l  0 ^ *C=I  j=-l  O & 


(3.23) 


where 


-Xr»;;L  _w  ced w: i ?+■»_ 

C©4A  rn<,  t\ 


- 


_ n-dod  HU  ({J-t-fc'^  4 (^T  Co,  h) 

C4<k3  »np 


(3.24) 
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Now  for  Equation  (3.14),  (3.20)  and  (3.23),  we  have  alto 

gether  NT  linear  algebraic  equations  for  NT  unknowns  which 

can  be  solved  simultaneously.  The  remaining  problem  we  face 
now  is  the  efficient  evaluation  of  the  integrals  defined 

by  in  (3.9),  (3.10)  and  (3.15).  the  P^'s  and  Q^'s  can  be 

evaluated  analytically  in  a rather  straight-forward  manner: 


% ( % ’7i ; V • V = I 1 ’ ** 

% 

tirfe-ljlfejkA-JjlVtJr'Tj)1)  t j Oj  -ty.) 


- V)V|  * -V>1  + Vi 


- “"S  - 

(3.25a) 


For  3*  and  simple  closed-form  solutions  appear  not 

possible.  While  (3.16a)  does  not  pose  too  much  of  a problem 
for  numerical  quadrature,  (3.16b)  have  rapidly  oscillatory 
integrands,  especially  when  m^  becomes  large.  We  found  that 


I 


Vi 


these  source  and  normal-dipole  integrals  having  oscillatory 
behavior  can  actually  be  written  in  terms  of  standard 
exponential  integral  with  complex  argument.  To  avoid 
disruption  of  the  continuity  of  exposition  of  the  main  text, 
we  provide  all  details  of  the  derivation  in  Appendix  A,  which 
is  based  on  a contour-interpretation  of  this  integrals  in 
the  complex  planet  The  results  can  be  summerized  as  follows. 

Define  the  complex  quantities: 

, ia.  = + <3’25> 

then 


1 'V—  t I MW. 

= -7T  e"  ^ 

* J ( e*1  e,(§2)  + € e'mf'e^'E1(|,)-eme5,£,e53)] 


(3.26) 


| 0 

Gru.(/'%.p  = ^7 

= .Ain  mi  - I ^ 


edjEI(3»)-e5*e,(-50-e  V,E,(5.)  + e’we' *'£,W 


(3.27) 


where  the  exponential  integral  function  E-^z)  is  defined 
to  be 

*The  author  is  grateful  to  Prof.  Yeung  for  providing  him  a 
copy  of  such  derivation. 

- 31  - 


I 


I 


where  ( = .5772156643  ) is  the  Euler's  constant  and 
double-precision  calculation  is  required  in  region  II  because 
of  slow  convergence.  For  region  III,  we  introduce  Laguerre 
quadrature  method  as  follows: 


e £,(^)  = L -<I: 


where 


i,  - (V-5? 


, I*=fef 2 ‘‘r 


It  is  well  known  that,  for  any  n, 


1 = e = z Ax  {(xT1)  = a 

where  the  xfn^  are  the  zeros  of  the  Laguerre  polynomial  Ln(t) 
and  A|n^  are  the  corresponding  Christoffel  numbers,  or 
weights.  Moreover,  it  is  known  (Szego,1939)  that  for  some 
the  error 


R = | I - Q | = (n  l) 

We  can  estimate  R by 

It? | < ("  ■ 1 L 

K ~ >o, 

1^1  - 2’rfl  2<  < o . 

Hence  one  sees  that  difficulties  in  using  the  Laguerre 
quadrature  scheme  are  encountered  when  z approaches  the 
negative  axis.  To  obtain  a small  value  of  |R|  , we  have 

to  make  an  optimum  choice  of  n,  which  is  about  |z|  in  the 
first  case,  |y|  in  the  second.  With  this  choice  of  n,  we 
find,  using  Stirling's  formula,  that  the  value  of  our  bounds 
for  i R I are  about  2'tTe  in  each  case. 
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VI.  SCATTERING  PROBLEMS 


VI . 1 Treatment  Of  The  Problem 

If  we  consider  the  scattering  problem  corresponding  to 
an  incoming  wave  system  from  the  left  with  unit  amplitude, 
the  definition  introduced  in  (2.1)  for  the  complex  function 
Q is  still  applicable  provided  one  sets  a^=0  and  a2=l/c* 
The  total  velocity  potential,  known  as  the  scattering 
potential,  is  normally  decomposed  into 

% = CPX  + <Vo  <4.1) 


where  CPX  , the  incident  wave  spatial  potential,  is  defined 
as : 


c p = 1A_  e+w 

1 t&aX,  nn0~-£' 


Aj  being  the  incident  wave  amplitude,  and  CPD  is  the  so- 
called  diffraction  potential  which  usually  has  to  be  solved 
for.  In  our  formulation,  it  is  much  more  convenient  to 
work  with  the  scattering  potential  directly.  The  scatter- 
ing potential  (4.1)  satisfies  the  boundary  conditions  (2.3b) 
and  (2.3c),  but  for  (2.3d) 

(<?$)„  = 0 on  S0  (4.3) 

So  far,  we  see  that  all  boundary  conditions  for  in 

Region  I are  homogeneous.  The  inhomogeneous  boundary 


i 
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condition  clearly  will  come  from  the  radiation  boundary. 

To  see  this,  we  observe  that  equation  (3.7)  is  no  longer  a 
complete  representation  of  the  solution  in  Region  II.  If 
the  right-moving  indident  wave  is  of  amplitude  A^  in  Region 


III,  (3.7)  should  be  replaced  by 


060  t j(f€~) 

t',  rru-ff 
(4.4) 


although  (3.6)  is  still  applicable,  viz.. 


(4.5, 


Now,  because  of  (4.3)  the  last  integral  of  (3.4) 
vanishes.  As  the  analog  of  (3.11),  we  obtain  the  following 
equation : 

-7T^P(f)+^  Cf  ^ r dA  + j l?o<j  r - y>j^  r)  4/5 

•^>L,Sb  i$F 


+ ZC<k  + t)  =-4e’ 


(4.6a) 


The  following  equations  are  the  results  when  the  field  point 
P is  placed  at  each  junction  boundary: 

[ +|  CP(^^r"V^r)^ 

I caP'ZA  rt\.  £ *=«  c^Wk-K. 


J 
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< I 


1 


tfi. 


^X/>-  to  in  i t/*  />•  h 


= -Aze  l(-Jo^>^+^^-3)]  for  p 


e Z, 


(4.6b) 


and 


l,cf>V,lJ  + L 

\Se> 


TCO»tyyC0l^' 

c*ol  yy£i~ 


<=i 


_ttcWk^H-)_  qT(0,m) 


£orpes- 


(4.6c) 


At  this  point  we  can  proceed  to  discretize  the  integral 
equation  as  we  have  previously  done  for  the  radiation 
problems.  Once  the  system  of  equations  is  sloved,  the 
complex-valued  reflection  R and  transmission  coefficient 
T are  simply  obtained  by 


R " A~/ax  * r " Ao/at 


(4.7) 


(4.8) 


The  phase  angles  are  defined  by 

arj  R = &R  = [~?m(9/Rc(R)) 

<u-jT  = 5"t  = ia»~'{J"(jyZe(T)] 

wt  ere  the  arc-tangent  is  understood  to  take  a range  of,i.e., 
-pj-  ^ .In  this  problem,  there  exist  two  sets  of 
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reflection  and  transmission  coefficients,  corresponding  to 
two  different  directions  of  incidence  of  the  propagating 
waves.  The  two  sets  are  somewhat  related  as  indicated  by 


Kreisel (1949)  and  Newman (1965) . If  we  denote  (R^,T^), 


(1*2^2)  as  the  reflection  and  transmission  coefficients  due 


to  the  incidence  of  the  unit-amplitude  propagating  wave  from 
the  left  and  the  right  respectively,  we  can  see  simply  write 
the  scattering  potential  as  follows; 
ml*  o 


9, 


(e^*+  R.e'*1"0*)*  (fl)  , * 


T,  e +cj) 


and 


<Pz 


( ( tB* m°  ^ 0 


+ 


R2e"m:/  ) 


l Tze^K  (p 


, % > - CO 


where 


Y±<  u\  = Ca^JL  mir  ( It-X") 

* ^ to**?? 


then 


(?,  = f?*  , 


RiR*  + I TT’CP/P')-  > 


(4.9) 


(4.10) 


(4.11) 


it  = _&*_  , RZR*  + T2T/(  1 

T«  R* 

0TZ  = P+  Ti 

where  rs*  — 4 — Z1W9  (4.12) 


I 
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and  the  asteriks  indicate  complex  conjucates.  These  are 
useful  relations  for  checking  the  consistency  of  the  numeri- 
cal solution. 


VI. 2 Relation  Between  Radiation  and  Scattering  Probelms 


So  far,  we  have  talked  about  the  'radiation'  and 


'scattering'  problems  seperately,  but  we  can  see  that  because 
of  the  mathematical  similarities  of  these  two  types  of 


probelms,  they  must  be  closely  related.  Certain  number  of 
interesting  identities  have  been  derived  in  a recent  work 


by  Newman (1975) , pertaining  to  the  case  of  equal  asymptotic 
depth.  Yeung (1975)  extended  results  of  these  identities  for 
the  case  of  unequal  depths.  The  details  are  reproduced  in 
Appendix  B for  completeness.  The  major  findings  are  sum- 
marized here.  If  A+,  B+  denote  the  asymptotic  wave- 
amplitudes  of  two  independent  solution  of  the  radiation 


problem,  then 


= i f aTb_ -ax 

. t;  ! A-d+  ~ A+  B-  (/}.B*  -A*  13. ) P/p-' , 


(4.13) 


which  can  be  used  to  obtain  the  coefficients  (R^,T^). 
(R2,T2)  then  follows  from  (4.11). 

Another  useful  finding  is  that  Haskind  relation  which 
relates  the  far  field  behavior  of  the  radiation  problem  to 
the  exciting  force  of  the  scattering  problem  can  be 
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generalized  to  the  case  of  unequal  depths  as  follows: 


F*  = -<:£VptYj±  (4.14) 

where  F*  is  the  exciting  force  or  moment  in  the  j-th 
direction  due  to  a unit-amplitude  incident  wave  and  Y±  is 
the  complex-wave  amplitude  as  x — y±<=o  corresponding  to  the 
forced  oscillation  of  the  body  or  any  part  of  the  bottom 
contour  in  the  j-th  direction. 


V.  RESULTS  AND  DISCUSSION 


The  method  discussed  in  section  III  and  IV,  has  been 
applied  to  a number  of  radiation  and  scattering  problems. 

In  all  cases  where  existing  results  are  available  and  are 
known  to  be  correct,  the  present  computation  shows  very 
good  agreement. 

V. 1 Radiation  Problems 

(i)  Circular  Cylinder  in  Heave  and  Sway 

The  problem  of  a circular  cylinder  in  heave  was 
first  solved  by  Ursell  (194 9) . It  has  become  a common  test 
problem  for  newly  develop  method.  Numerical  solution  for 
a heaving  and  sway  circular  cylinder  have  been  obrained  for 

frequencies,  K=  -?-A  = 0.  to  2.0,  with  cl  being  the  radius  of 

a 

the  circle,  for  purpose  of  comparision  with  others.  Table  1 
shows  the  sway  added-mass  and  damping  coefficient  (as  well 
as  the  wave-amplitude  ratio)  compared  with  Potash (1970)  who 
used  the  Havelock  wave  source  to  solve  the  problem.  To 
simulate  the  infinite-depth  situation,  we  exploit  the  fact 
that  depth  effects  for  any  wave  decreases  approximately  in 
the  form  of  e m°\  For  a given  frequency,  we  have  chosen 
h =^A  • The  agreement  is  clearly  assertive  of  the 

correctness  of  the  programming  efforts. 

Results  for  the  heave  added-mass  and  damping  coeffi- 
cients in  water  of  finite  but  constant  depths  are  shown  in 
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r; 

r"  ~ 

Kim 

Potash 

1.147 

1.156 

1.282 

1.290 

1.167 

1.174 

0.988 

0.993 

0.810 

0.812 

0.650 

0.659 

0.452 

0.450 

0.335 

0.332 

0.269 

0.265 

0.231 

0.226 

0.211 

0.205 

0.199 

0.193 

0.193 

0.187 

0.190 

0.185 

0.189 

0.186 

0.064 

0.499 

0.716 

0.843 

0.887 

0.877 

0.791 

0.688 

0.597 

0.520 

0.457 

0.405 

0.362 

0.324 

0.291 


0.063 

0.501 

0.712 

0.850 

0.895 

0.885 

0.798 

0.695 

0.604 

0.526 

0.462 

0.408 

0.363 

0.324 

0.291 


a Wave  Amp 

'"'ll-  ^ 

M 

»+.on  Amp 

Kim 

Potash 

0.031 

0.032 

0.265 

0.266 

0.423 

0.426 

0.574 

0.578 

0.707 

0.711 

0.821 

0.825 

1.003 

1.008 

1.144 

1.150 

1.259 

1.266 

1.356 

1.364 

1.442 

1.448 

1.518 

1.521 

1.585 

1.585 

1.645 

1.641 

1.695 

1.691 

r of  points  on  body=19 


'able  1.  Added- 
Wave-A 


lg  Coefficients  and 
:io  in  Sway  Motion 


npl  it u d e 


Figure  3 and  compared  with  Bai  and  Yeung (1974).  In  this 
case,  h/a=2.0.  The  present  calculations  yield  a value  of 
heave  added-mass  coefficient  equal  to  0.500  at  v’d^O.OOl 
which  is  consistent  with  the  value  of  0.4984  obtained  by  an 
alternate  approach  of  considering  the  limiting  boundary- 
vlaue  problem  (Yeuncr  and  Hewman,  1976) 

Another  consistency  check  consists  of  the 
Equation (2 . 6)  relating  the  local  solution  obtained  by  inte- 
grating the  potential  on  the  body  and  the  far-field  wave 
amplitudes  defined  by  the  first  term  of  the  eigenfunction 
expansion  Eq.(3.6)  and  (3.7).  For  a fluid  of  unequal  depths 
on  each  side,  this  was  shown  by  Yeung  (1973)  to  be 

— D + 1 D'  = 1 (5.1) 

a o- 

/ 

This  equality  is  satisfied  within  0.5%  accuracy  in  our 
calculations . 

As  a gage  on  the  convergence  characteristics  of  the 

solution  as  a function  of  the  location  of  the  radiation  boundary, 
or  the  number  of  terms  used  in  the  eigenfunction  expan- 

sions, a small  numerical  experiment  was  created  and  results 
are  summarized  in  Table  2 for  va=0. 9 and  h/a=5.0. 

The  final  algorithm  developed  for  the  choice  of  %K  or 
versus  N„  or  NT  is  as  follows: 
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Number  of  Eigenfunctions  = 5 


5r 

A 2 z 

X 

2.745 

0.5937 

0.4491 

0.7561 

3.000 

0.5959 

0.4490 

0.7561 

4.500 

0.5991 

0.4582 

0.7642 

8.000 

0.5996 

0.4579 

0.7640 

15.000 

0.5997 

0.4578 

0.7639 

Location  of  Junction  Boundary ( ) = 3.0 


Nr 

— 

1 

0.6448 

0.3589 



0.6414 

3 

0.6081 

0.4686 

0.7780 

5 

0.5959 

0.4490 

0.7561 

10 

0.5961 

0.4451 

— 

0.7552 

Table  2.  Summary  of  Results  for  Varying  the  Location  of 

Junction  Boundary  or  the  Number  of  Eigenfunctions 

We  consider  two  different  cases,  one  is  the  uniform  depth  and 
the  other  the  varying  depth  case.  For  the  uniform  depth 
case,  in  general,  taking  the  junction  boundary  or 

as  close  to  the  body  as  possible  will  reduce  the  computa- 
tional time  substantially.  But  in  this  case,  local  distur- 
bances play  such  an  important  role  that  a large  number  of 
eigenfunctions  are  required  to  treat  the  problem  properly. 

To  avoid  this  deficiency,  we  usually  specify  about  one-third  of 
a wave-length  from  the  body  as  an  input  value  of  or  for 
each  frequency,  and  2J)  is  considered  as  the  maximum  number  of 
eigenfunctions.  As  we  mentioned  in  Section  III. 2,  the  number 
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of  terms  can  be  less  than  20  if  the  decay  factor  d(K,?)=<2 
satisfies  the  error  criteria  before  ND  or  N has  the  value  20. 
If  not,  we  can  change  the  distance  of  the  junction  boundary 
from  the  body  until  the  decay  factor  satisfies  the  error 
criteria  with  NR  or  set  equal  to  20.  Concerning  the  vary- 
ing depth  case,  we  can  use  exactly  the  same  scheme  as 
the  above  case  except  for  the  initial  choice  of  or  . 

(ii)  Rectangular  Cylinder  in  Heave  and  Sway 

The  heave  and  sway  added-mass  and  damping  coeffi- 
cients in  water  of  finite  but  constant  depth  are  shown  in 
Figure  4.  These  results  are  compared  with  Lebreton  and 
Margnac  (1966)  who  used  the  Green- function  integral  equation 
technique  to  solve  the  problem.  The  limiting  added-mass 
coefficient  is  0.822  as  the  frequency  goes  to  zero. 

(iii)  Bulbous  Section  in  Heave  or  Sway 

The  bulbous  section  being  considered  consists  of 
a fully  submerged  circle  with  a narrow  vertical  stem  mounted 
on  top  of  it  as  shown  in  Figure  5.  Here,  we  considered 
various  cases  of  depth  changed,  h+=  l.ld,  1.25d,  1.50d  and 
2 . Od  holding  h = 2.0d.  Results  for  the  force  coefficients 
due  to  sway  motion  are  given  in  Figure  6. a,  and  those  due 
to  heave  motion,  in  Figure  6.b.  It  is  of  interest  to  note 
that  for  the  heave  motion,  the  added-mass  coefficients 
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Fig.  5 Bottom  Geometry  with  Bulbous  Section 

increase  proportionally  to  a decrease  of  depth  throughout 
most  of  the  frequency  range,  but  the  damping  coefficients 
remain  practically  the  same,  especially  when  va  0.15. 

On  the  other  hand,  hydrodynamic  coefficients  associated  with 
sway  mode  behave  differently.  In  all  cases,  when  va  fc  0.20 
coefficients  are  almost  identical.  One  phenomenon  we  notice 
is  that  in  some  frequency  range  the  sway  added-masses  become 
vanishingly  small  and  even  slightly  negative.  Such  an  occu- 
rance  of  negative  sway  added-mass  observed  by  Ogilvie (1963 ) 
when  a fully  submerged  circular  cylinder  undergoes  sway. 
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Fig  6a  Added  - Mass  and  Damping  Coefficient  for  Bulbous  Section  in  Sway 


V. 2 Scattering  Problems 
(i)  Sinusoidal  Hump 


The  geometry  and  results  are  shown  in  Figure  7. 

Yt(x)  is  the  total  wave-amplitude  function  due  to  an  incident 
wave  of  unit  amplitude  propagating  from  the  left  over  the 
sinusoidal  hump.  The  results  obtained  here  show  very  good 
agreement  with  Bai(1972)  obtained  using  a finite  element 
method.  The  complex  wave-function  YT(x)  is  defined  as 

= e YpW  * Y,  (70  + * X C<)  (5.2) 

where  the  first  term  after  the  first  equal  sign  represents 
the  incident  wave  with  unit-amplitude  and  the  second  the 
diffracted  wave  profile.  This  test  essentially  verifies 
our  theory  discussed  in  Section  IV. 1 for  the  treatment  of 
the  scattering  problem. 

(ii)  Finite-Size  Step 

The  reflection  and  transmission  coefficients  of 
a finite-size  step  are  given  in  Table  3 and  4,  and  compared 
with  Hilaly(1967)  for  h /h+  = 2.5  and  6.5.  From  equation 
(4.8)  and  (4.11),  we  expect 

|R,l  = |Rzl  (5-3) 

I R.r+  nf  ( p+p~)  = '•  <5-4> 

where  (R^,T^)  and  (R2>T2)  are  the  resulting  reflection  and 
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Eqn  (5.4) 


vh 

R 

1 

T1 

Eqn  (5.4' 

!Ril 

arg  Rx  (°  ) 

lTil 

arg  Tj_  (°  ) 

L . H . S 

0.03395 

0.4269 

-163.14 

0.5799 

4.99 

1.00025 

0.4273* 

-163.26* 

0.5797* 

4.97* 

1.00001 

0.04365 

0.4242 

-160.86 

0.5846 

5.64 

1.00016 

0.4246 

-161.00 

0.5846 

5.61 

1.00001 

0.05813 

0.4201 

-157.86 

0.5918 

6.46 

1.00011 

0.4205 

-158.02 

0.5917 

6.44 

1.00005 

0.06644 

0.4134 

-153.73 

0.6035 

7.56 

1.00002 

0.4138 

-153.93 

0.6034 

7.53 

1.00007 

0.12115 

0.4016 

-147.68 

0.6244 

9.04 

0.99980 

0.4021 

-147.94 

0.6244 

9.02 

1.00008 

0.20031 

0.3774 

-137.97 

0.6679 

11.07 

0.99938 

0.39260 

0.3158 

-120.25 

0.7803 

13.05 

0.99851 

0.3162 

-120.83 

0.7810 

13.03 

1.00023 

1.09093 

0.1458 

-90.67 

1.0395 

7.61 

0.99760 

0.1470 

-91.99 

1.0407 

7.69 

1.00025 

* Taken  from  Hilaly (1967)  . His  phase  angle  is  related  to 
the  present  notation  by:  a = - ( tt+  arg  R) 


UNIT 


Table  3.  Reflection  and  Transmission 
Coefficients  of  a Step 


CASE 

I 

R1 

arg  Rj 

T1 

arg  Tx 

2212 

-170.6 

.7820 

1.65 

223* 

-170.7* 

. 7802* 

1.66* 

2206 

-169.3 

. 7835 

1.86 

2223 

-169.4 

.7818 

1.88 

2196 

-167.6 

. 7858 

2.14 

2213 

-167.8 

. 7841 

2.17 

2181 

-165.3 

. 7895 

2.52 

2198 

-165.5 

.7878 

2.56 

2154 

-162.0 

. 7960 

3.05 

2170 

-162.3 

. 7944 

3.01 

2100 

-156.6 

. 8091 

3.85 

2114 

-157.0 

. 8076 

3.91 

1964 

-146.5 

. 8421 

5.08 

1973 

-147.1 

. 8410 

5.18 

1423 

-120.0 

. 9657 

5.91 

1413 

-121.1 

. 9654 

6.07 

CASE  II 
R2  arg  R2 


arg  T2 


* Taken  from  Hilaly.  His  phase  angle  is  related  to  the 
present  notation  by:  a = - (ir+arg  R) 


h+  h + 


— 2.5 


(Case  I) 


(Case  II) 


Table  4.  Comparision  of  Reflection  and  Transmission 
Coefficients  for  Incidents  Waves  in 
Opposite  Directions 


transmission  coefficients  when  an  incident  wave  of  unit- 


amplitude  comes  from  the  left  and  right  respectively.  The 
ratio  D+/D  evidently  represents  the  group-velocity  ratio. 

The  phase  angles  given  in  these  tables  are  defined  by  (4.8). 
Our  results  not  only  show  very  good  agreement  with  Hilaly's, 
but  also  meet  the  testing  of  Equation  (5.3)  and  (5.4)  with 
sufficient  accuracy. 

(iii)  Steep  Hump  at  a Step 

We  apply  our  numerical  method  to  determine  the 
reflection  ability  of  a steep  hump  located  in  front  of  a 
step.  This  is  one  sample  of  physical  geometry  which  can 
not  be  easily  accomodated  by  any  analytical  methods.  The 
geometry  of  the  hump  and  the  results  are  presented  in  Figure 
8.  From  Equation  (4.11)  one  can  define  an  energy  trans- 
mission  factor  k by  fc  = | T I D /D  . The  hydrodynamic 
characteristics,  reflection  and  transmission  coefficients, 
associated  with  a finite-size  step  without  the  hump (Table  4) 
are  also  plotted  for  the  purpose  of  comparision.  One  notices 
that  the  hump  causes  a considerable  increase  in  wave 
reflection,  and  the  energy  transmission  factor  achieves  a 
minimum  at  a frequency  of  vh  =1.0. 


Transmission  Coeffs. 


1 


VI . REMARKS 


In  this  thesis,  we  have  presented  a versatile  numerical 
method  of  solution  which  is  applicable  to  the  two-dimensional 
time-harmonic  free-surface  flow  problem  for  a fluid  of  arbi- 
trary bottom  topography  and  have  investigated  the  validity  of 
this  method  by  comparing  our  results  with  others  known  to  be 
correct.  In  every  case,  we  obtained  very  satisfactory 
agreement.  One  item  we  should  mention  is  that  for  scattering 
problems  the  numerical  accuracy  is  not  as  good  as  expected 
when  the  reflection  coefficient  becomes  vanishingly  small  as 
the  frequency  becomes  fairly  high.  To  avoid  this  deficiency, 
it  seems  to  be  worthwhile  to  make  a modification  in  computa- 
tional procedure.  It  is  possible  to  extend  this  method  to 
a general  three-dimensional  problem  without  any  appreciable 
amount  of  conceptual  difficulties. 
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APPENDIX  A 


The  Integrals  J'k  and  Gk 


those  defined  by  (3.10): 


l° , J-whr 

= Jrj  COd  0 

1 w V 


i 


We  will  show  that  they  can  be  written  in  terms  of  exponential 
integrals.  The  derivation  to  follow  is  considerably  simpler 
if  we  introduce  the  following  complex  variables: 

^ = oM  ^ , oL  = ^ 

$ = T + *1 

where  for  practical  purposes  one  only  has  to  consider  the  case 
(;x-$)<o  since  Cr  remains  the  same  while  J differs  only  by 
a sign  for  the  opposite  case.  In  the  complex  plane  Gr(z)  is 
given  by: 

G (d-g)  = f?e  \c  (W.  m (£+<£)  w (i-*)) 

= (A. 2) 

where  the  contour  of  integration  CQ  is  a vertical  line  from 
% = to  <^=o  as  shown  in  Fig.  A.l. 


i 


If  W ( z ) is  known,  the  evaluation  of  3 is  straightforward 


since 


which  follows  from  the  fact  that  the  derivative  of  any  analy- 
tic function  is  independent  of  the  direction  of  approach. 

Now  for  W(z),  an  integration  by  parts  yields 


im  ( K + 


~<5-Llvl  (X.m'Q.)  hi  ^ + 2,  C 


l ( 


da 


TnU+A-fc) 


+ 


-rr\(-^+^£) 


- wt^+.C'iv) 


(A. 4) 


The  paths  of  the  integrals  of  (A. 4)  in  the  u-plane  are  shown 


in  Fig.  A. 2. 


\jw  Cii) 


K-mi 

\ 

\ 


Ke(u) 
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Both  integrals  can  be  replaced  by  a pair  of  integrals  exten- 


. 

J 


! 


1 


ding  from  the  end  points  to  z=  +<?°  . Exploiting  the  defini- 


tion of  the  exponential  integral  E^(z) 


-t 

I"PI 

IT 


(A. 5) 


where  the  Riemann  cut  is  along  the  negative  real  axis,  we 
obtain  for  (A.4),  after  accounting  for  the  pole  at  u=0  for 
the  first  integral,  the  following  expression: 


Wc-p  - < I ( iml)  + e.  ^ ^-E.  t E,  ( m(}+.<;£))  +2tt) 


I v | 


Whence , 

GfC®^)  = -dX*  ml  U U'+p*  - V c^s  w>(Jj+£) 


(A. 6) 
Wo( 


- e"<^+‘'f,)[El(-mJ)-E,(-m(j+^W]  j (A  7) 


which  is  (3.27).  Next  for  we  recall  that  by  Liebnitz 

rule 


T^Eitp  - eie,(i)~  f 


(A. 8) 


Using  (A. 3)  and  (A. 6),  we  have  for  the  normal-derivative 
integral : 


- m (oi) 


3(d.Jp  = -TT^jn(oi)  e ccKiniyl) 

tjl  E,  (vnly^)  - E,(w$j) 
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r*." ' 


+ gW(i+^-R)re, El(_w^))i 

(A. 9) 

which  is  identical  to  that  indicated  in  (3.26).  Note  that 
as  o(->o  the  exponential-integral  terms  vanish  in  pairs 
because 

r* 


(A. 10) 


and  (A. 9)  reduces  to  merely 
± 


$ ( 6 * 3)  = (°0  (a. id 

which  follows  naturally  from  the  physical  interpretation  of 
the  integral  as  the  normal  velocity  of  a source— distribution . 
In  this  same  limit,  C^W.H)  can  be  written  as  the  sine-  and 
cosine-integrals : 

6rC0.Jp  = Awi  nry^  JPoj  1^1  t Axm  (WJ+fr)  ~Ci  (“"»})] 


(A. 12) 

where  C.  and  are  those  defined  in  Abramowitz  and  Stegun 
(1970) . 
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APPENDIX  B 

Symmetry  Relations  for  Radiation  and  Scattering  Problems 

Let  ^ and  ^ be  the  spatial  scattering  potentials  co- 
rresponding to  a unit  amplitude  right-  or  left-moving  incident 
wave  respectively.  Therefore  asymptotically 


<P, 


i C e 


T,  ex'<X  * + 1 


(B.l) 


and 


(B.  2 ) 


Klp>  / 


where  the  R’s  and  T’s  are  complex-valued  transmission  and 
reflection  coefficients. 

Next,  define  the  integral  operator  on  two  poten- 

tial functions  cp  and  as 


\ 


<P.W>  = j l 


(B . 3 ) 


where  C consists  of  the  boundary  contours  S_f  S_,,  S . and 

r Jd  O 

two  vertical  boundaries  at  infinities.  Making  use  of  the 
boundary  conditions  (2.3b),  (2.3c),  and  (4.3)  as  well  as  the 

asymptotic  behavior  of  and  Cf> * we  have 


<Cp,  , cp,*>  : 1 - p,  R*  t ( T.T* ) D%T 


(B . 4 ) 
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< q>2 , 0 : 1 = f?2  Ra*  + (T*  t,*)  v/d* 


<<P. , <P*>:  t2  D~  = ID1 


(B.6) 


cp, , <&*>:  Tz  R,  0~  - T,  R/p- 


(B.  7) 


where  D*  (m“h4  ) is  defined  earlier  by  (2.18).  The  first  two 
equations  are  statements  of  conservation  of  energy.  Together 
with  (B.6),  one  can  deduce  easily  that 


IP, I = IPzl 


(B.  8 ) 


(B. 7 ) and  (B.8)  can  also  be  combined  to  yield 


T,  _ R» 
T R. 


(B . 9 ) 


Hence  if  T^  and  R^  are  known,  (B.6),  (B.8),  and  (B.9)  define 

R2  and  T2  completely.  These  are  essentially  restatements 
of  observation  by  Kreisel  (1949)  and  Newman (1965) . We  reproduce 
them  here  for  the  sake  of  completeness. 

To  obtain  relations  between  the  radiation  problems  and 
scattering  probelms,  we  denote  the  radiation  potential  by  , 
which  has  an  asymptotic  behavior  of 


ax<t 


A_  X"(J>  e 


7C  — v -+ 

7<  • — *■  - <30 


(B. 10) 


and  boundary  condition  on  Sq  defined  by  (2.3d).  Now  applying 
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the  operator  (B.3)  to  the  following  combination  of  and  CfJ , 
or  (-Pz  we  have 


<<P,X> : 
<A  q>>: 


D/X  t d/4*T,  + ^ = o 

r 

D A -X  CJ>-fu)eU  = o 


Sa 

DAT/  + DX  R* 


<cpz,op>:  DX  <?2{m)cU  = O 


(B. 11) 
(B. 12 ) 
(B. 13) 
(B . 14 ) 


This  basically  covers  all  possible  combinations.  Note  that 
(B.12)  and  (B.14)  are  essentially  Haskind  relations  for  the 
unequal  depth  case;  namely  if  denotes  the  excitation  force 
or  moment  on  the  body  Sq  in  the  j-th  direction  due  to  unit- 
amplitude  incident  wave  then 


p/  = D*(nifc)Y0*  (B.  15) 

A 

where  Y*  are  the  asymptotic  wave  amplitude  of  cp  with  a 
boundary  conditions  f(s)  corresponding  to  a unit  oscillation 
of  the  j-th  mode.  (B.ll)  and  (B.12)  can  be  combined  together 
noting  that  f (s)  is  real  for  radiation  problems;  whence 

DA*  R,  + D/CT i + D~A-  ~ 0 (B.16) 

Similary,  for  (B.13)  and  (B.14) 

DAT/  + D XR/  + D*A*  = o (b.  17 ) 
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It  is  not  obvious  that  (B.17)  is  a redundant  relation.  To 

★ 

show  this,  we  multiply  (B.16)  by  and  add  the  result  to 
itself, 

DA*  R.a  + D/Ct,  r:  ~ D+  A+T,  * - DA*  = o (B.  is) 


Next,  making  use  of  (B.4)  for  the  first  term  we  have 


( - D+/C  T,  + D+/i 


~n  r, 
+ T,# 


DA.  T,  = ° 


(B . 19 ) 


Hence  by  (B.6)  and  (B.9),  we  recover  (B.16).  Therefore,  as 
in  the  constant  depth  case,  only  one  equation  is  available. 
But  unlike  the  case  considered  by  Newman  (1 975 ) , because  of 
the  bottom  geometry,  a symmetrical  body  will  not  necessarily 
generate  symmetric  or  anti-symmetrical  wave  motion.  Similar 
simplications  do  not  result.  On  the  other  hand,  deduction 
of  the  value  of  transmission  and  reflection  coefficients  from 
the  radiation  problem  is  still  possible. 

Denote  the  modified  wave  amplitude  at  infinities  of 
two  modes  of  oscillations  by  /4+,  B+ » then 


At  - D Ya 
B+  " D Y» 


(B. 20) 


where  Y^,  Y*  are  the  complex  spatial  wave  amplitude  defined 
by  (2.8d)  with  |a|  =1.0.  Straightforward  application  of 
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